Skip to content

Astra curved detector - #1729

Open
NogginBops wants to merge 7 commits into
odlgroup:masterfrom
NogginBops:astra-curved-detector
Open

Astra curved detector#1729
NogginBops wants to merge 7 commits into
odlgroup:masterfrom
NogginBops:astra-curved-detector

Conversation

@NogginBops

@NogginBops NogginBops commented May 5, 2026

Copy link
Copy Markdown

This PR is based on #1634 but remade for odl 1.0.0.

This PR uses cyl_cone_vec from ASTRA 2.4.0 to allow curved cone beam detectors.

From #1634 there was a note about a off-by-half bug related to pixel centers that I'm not sure if it still applies.

Known bug: there is a misinterpretion of the detector center if the detector partition isn't centered around 0, which is the case for quarter-pixel-shifted detectors. Fixing this needs to be coordinated with the Mayo CT data loader, which currently is assuming the presence of this bug.

There is also a FIXME in astra_cuda.py related to not being able to arbitrarily transposing pytorch tensors. Not sure what the best solution for that is, guidance is appreciated. :) This is fixed by using permute_dims.

@leftaroundabout could you take a look at and see if you can see anything that is obviously wrong that I could fix? You can do a full review when you have time.

Comment thread src/odl/applications/tomo/backends/astra_cuda.py Outdated
@NogginBops

Copy link
Copy Markdown
Author

There is a bug in this code that I will investigate tomorrow, but because creating pull request threads is down (https://www.githubstatus.com/) can't comment this on the affected line of code.

On line 253 in astra_cuda.py:
I have observed a crash here that I will investigate tomorrow:

Traceback (most recent call last):
  File "/workspace/reconstruction.py", line 182, in <module>
    main()
  File "/workspace/reconstruction.py", line 133, in main
    resid = sinogram - A(x)   # residual in data space
                       ^^^^
  File "/opt/odl/src/odl/core/operator/operator.py", line 663, in __call__
    out = self._call_out_of_place(x, **kwargs)
          ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/opt/odl/src/odl/applications/tomo/operators/ray_trafo.py", line 304, in _call
    return self.get_impl(self.use_cache).call_forward(x, out, **kwargs)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/opt/odl/src/odl/applications/tomo/backends/util.py", line 47, in wrapper
    return fn(self, x, out, **kwargs)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/opt/odl/src/odl/applications/tomo/backends/astra_cuda.py", line 175, in call_forward
    return self._call_forward_real(x, out, **kwargs)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/opt/odl/src/odl/applications/tomo/backends/astra_cuda.py", line 253, in _call_forward_real
    return self.proj_space.element(proj_data)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/opt/odl/src/odl/core/discr/discr_space.py", line 369, in element
    self, self.tspace.element(inp, copy=copy)
          ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/opt/odl/src/odl/core/space/base_tensors.py", line 686, in element
    return wrapped_array(arr)
           ^^^^^^^^^^^^^^^^^^
  File "/opt/odl/src/odl/core/space/base_tensors.py", line 625, in wrapped_array
    raise ValueError(
ValueError: shape of `inp` not equal to space shape: (1560, 230, 11066) != (11066, 1560, 230)

@leftaroundabout

Copy link
Copy Markdown
Contributor

@NogginBops looks like permute_dims behaves in a subtly different way to the old NumPy transpose. Not very surprising... in my experience this kind of thing always requires a lot of trial and error to get right in Python.
Good work regardless, getting it this far!

@NogginBops

Copy link
Copy Markdown
Author

I've seemingly fixed the issue. The issue was that the old permutation had order 2 which meant applying the same permutation again returned to the original shape, but the new permutation is order 3 which means applying the same permutation again does not return to the dimentions to the original order. The fix is to introduce explicitly specify the inverse permutation.

@leftaroundabout

Copy link
Copy Markdown
Contributor

@NogginBops sounds reasonable. Thanks a lot for your efforts!

I don't see anything wrong with the code, but still need to do a bunch of testing before I can merge it.

@leftaroundabout

Copy link
Copy Markdown
Contributor

Well, this backprojection seems to suggest that at least the axes are correctly mapped:

image
detector_partition = odl.uniform_partition([-np.pi/8, -160], [np.pi/8, 160], [512, 512])
geometry = odl.applications.tomo.ConeBeamGeometry(
    angle_partition, detector_partition, src_radius=200,
    det_radius=200, det_curvature_radius=(400, None),
    axis=[1, 0, 0])

(everything else the same as in ray_trafo_cone_3d).

What's strange is that I have to make the detector very "oversized" (so the projection of the actual phantom only takes up a small part of the data space, which would be rather inefficient). If not, then the backprojection comes out heavily cropped on the sides, even though the projection is still well within the detector area:

image

Not sure if this is because of some actual bug in the code, or I'm just interpreting the results wrong.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants